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>■ ■ Abstract 

o ■ 

I We discuss two approaches to the calculation of matrix elements of the Argonne vis 

• potential. The first approach is applicable in the case of a single-particle basis of 

harmonic-oscillator wave functions. In this case we use the Talmi transformation, 
implemented numerically using the Moshinsky transformation brackets, to separate 
the center-of-mass and relative coordinates degrees of freedom. Integrals involving 
the radial part of the potential are performed using Gauss-Hermite quadrature 
y ' formulas, and convergence is achieved for sets of at least 512 points. We validate 

' the calculation of matrix elements of the Argonne fig potential using a second 

approach suitable for the case of an arbitrary functional form of the single-particle 
wave functions. When the model space is represented in terms of harmonic-oscillator 
' wave functions, results obtained using these two approaches are shown to to be 

identical within numerical accuracy. 
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' 1 Introduction 



Direct comparison of experimental data and theory become ambiguous if prob- 
lems such as the many-body problem are not solved accurately. This is par- 
ticularly true in nuclear physics, where one uses numerical solutions of the 
Schrodinger equation to constrain in part the effective interaction between nu- 
cleons [1]. Despite tremendous accomplishments in the past decade in carrying 
out ab initio studies of nuclear structure of light systems [2|3P|5|6|7f8|9l[T0|llfl2fl3fl4fl5l 
numerical solutions of medium and large nuclear systems are yet to be per- 
formed in the framework of realistic nuclear interactions [TB1I17II18I. 
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The matrix elements of the interaction represent the building blocks of any 
many-body approach which shares features in common with the shell model 
approach. Here, the model space is represented in terms of a complete set 
of single-particle wave functions. All possible single-particle configurations 
with the appropriate properties should in principle be used to diagonalize the 
Hamiltonian, and obtain the associated spectrum of eigenvalues and eigenfunc- 
tions. However, the single-particle model space is in principle infinite. Hence, 
finite reductions of the model space must be considered and convergence to the 
continuum limit must be achieved. The efficiency of this procedure is linked 
to the ability of the functional form of the single-particle wave functions to 
capture the relevant correlations at both short- and long-range scale. 

Traditionally, matrix elements of the interaction are formed in a basis of 
harmonic-oscillator wave functions. This has the advantage that one can per- 
form an exact separation of the center-of-mass and relative coordinates degrees 
of freedom using the Talmi transformation approach [12] (via its numerical im- 
plementation based on the Moshinsky transformation brackets [2D]), and inte- 
grals can be accurately performed using Gauss-Hermite quadrature formulas. 
Therefore, we will by discussing the calculation of the two-body matrix ele- 
ments in a harmonic-oscillator basis. This is done both for completeness, and 
for testing purposes, as the more elaborated procedures to be discussed in the 
second part of this paper will require validation. 

Unfortunately, the asymptotic form of the harmonic-oscillator wave functions 
is Gaussian which leads to densities decaying as , whereas the natural 
tail of the nuclear density is exponential, e~^. Added flexibility is desirable, 
and can be achieved by considering linear combinations of harmonic oscilla- 
tor wave functions [T0llllfl2j . or different functional representations for the 
hole and particle sides of the spectrum [2T]. To this effect, in the second part 
of this paper, we will discuss the framework of calculating matrix elements 
of the Argonne vis potential in a basis of single-particle wave functions of 
arbitrary functional form. Here, the ability of using the Moshinsky transfor- 
mation brackets approach to perform the separation of the center-of-mass and 
relative coordinates degrees of freedom is rendered numerically unfeasible, as 
the required CPU time is of order compared to the calculation involving 
Moshinsky transformation brackets. 

For illustrative and testing purposes, we will apply this approach to the case 
when the radial part of the single-particle wave functions, Rni{x) can be ex- 
panded out into harmonic-oscillator wave functions, HOkiix), i.e. 



This is convenient because we can still use a Gaussian quadrature technique to 
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perform the numerical computation of the radial integrals. This also allows us 
to check of the reliability of our matrix-elements calculation by comparing with 
matrix elements calculated in a harmonic-oscillator basis using the procedures 
using the Moshinsky transformation brackets. 

The paper is organized as follows: In Sec. [2l we review features of the Ar- 
gonne wig potential, together with the definitions of the matrix elements to be 
discussed further. Basic aspects of the separation of the center-of-mass and 
relative coordinates using the numerical implementation of the Talmi trans- 
formation via the Moshinsky transformation-brackets approach, together with 
the procedures required for the 7 types of operators part of the Argonne vis 
potential, are discussed in Sec. |3l In Sec. |4] we present the general results 
used to calculate the matrix elements of a translationally-invariant two-body 
interaction in an arbitrary single-particle basis. The implementation of this 
general approach for the particular case of the 7 operators in the Argonne vis 
potential is discussed in Sec. [51 This calculation is carried out in a jj-coupling 
scheme. The part of the calculation involving the isospin part of the interac- 
tion is identical in both frameworks. A brief review of this calculation is the 
object of Sec. O Finally, in Sec. [7] we discuss numerical convergence aspects 
of the matrix elements calculation corresponding to a ^^0-like model space of 
harmonic-oscillator single-particle wave functions, and compare with results 
obtained using the two frameworks. We also include an appendix on the cal- 
culation of large sets of abscissas and weights for the evaluation of integrals 
using Gauss-Hermite quadrature formulas, followed by an appendix summa- 
rizing results relevant to the calculation of reduced matrix elements used in 



2 Preliminaries 

Consider the two-body approximation of the Hamiltonian describing a system 
consisting of protons and neutrons 



SecO 



(2) 



i i<j 



In the second- quantization representation this becomes 



H=E 



(3) 



where Greek letters label the single-particle states \a) = \nlsjmj; ^rriT-), with 
n > 1, s = |, j = / and = +|(— |) - for a proton (neutron). The parity 
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of these states is (— 1)^ In this work we will assume that the radial part of 
the single-particle wave functions, Rni{x), is simply given by the radial part of 
the harmonic oscillator wave functions T-LOkiix), subject to the normalization 
conditions 

/ [x'^ dx]Rnl{x) Rn'l'ix) = 6nn'Su' . (4) 

Jo 

Here, x = r/b, with b the oscillator parameter. We note that the single-particle 
radial functions 7^„i(r) are defined in terms of the the radial functions Rni{x) 
such that they satisfy the normalization condition 



oo 
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Thus, we have 



[r dr]7^„/(r)7^„/^/(r) = Snn'Sw ■ (5) 



= ^^r.l{x). (6) 



2.1 Argonne Vis potential 

The Argonne Vis potential [TB] is an updated version of the nonrelativistic 
Argonne potential [22] that fits both np data and pp data, as well as low- 
energy nn data scattering parameters and deuteron properties. The potential 
was fit directly to the Nijmegen A^A'^ scattering database, which contains 1787 
pp and 2514 np data in the range 0-350 MeV, together with the nn scattering 
length measured in d{7r~,'y)nn experiments and the deuteron binding energy. 
The fit has a P^r datum of 1.09. 

The strong interaction part of the Argonne vis potential is projected into an 
operator format with 18 terms: 

i) a charge-independent part that has 14 operator components (as in the older 
Argonne Ui4 [22]), featuring 7 basic operators, 

1, ai ■ aj, i ■ S, £\ f (a, ■ a,), {i ■ S)^ , (7) 
and their charge-independent counterparts 

n-Tj, {ai-(rj){n-Tj), Sij{Ti-rj), {i ■ S) ■ Tj) , 

e (r, ■ r,), e (a, ■ a,) (r, ■ r,), (£ ■ Sf (r, ■ r,) ; (8) 



ii) a charge-independence breaking part that has three charge-dependent op- 
erators 
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Tij, {(Ti ■ (Tj)Tij, Sij Tij , (9) 

where Tij = 3 T^iTzj — Ti ■ Tj is the isotensor operator, defined analogous to the 
Sij operator, and iii) one charge-asymmetric operator 



r,i + T,j . (10) 

In principle, there could be more charge-independence breaking terms, such as 
{i ■ S) Tij or Sij {jzi + Tzj), but the scattering data are not sufficiently precise 
to identify them at present. The potential also includes a complete electro- 
magnetic potential, containing Coulomb, Darwin-Foldy, vacuum polarization, 
and magnetic moment terms with finite-size effects. 



2.2 Two-body matrix elements 



Matrix elements of the two-body interaction can be specified either in particle- 
particle (pp) coupling or in particle-hole [ph) coupling. Either set of matrix 
elements completely specify the two-body interaction. We define ph-conpled 
matrix elements as 

m,\vm,)= y: (mi^^) (n) 

mim2m3m4 

X (_)i3-m3 {j^rnijs -mg | A/i) (-^^^^^ ^j^m^j^ -ms | A//) . 
Correspondingly, the j?p-coupled matrix elements are defined as 

((12)^1^1(34)^) = (12|\/|34)(jimu2m2 | L M){3^m^ j^m^ \ L M) . 

(12) 

In practice, the calculation of the two-body matrix elements using a harmonic- 
oscillator single-particle basis set is carried out in pp coupling. Subsequently, 
j»/i-coupled matrix elements are evaluated from their pp counterpart, using 
the relationship 

((13),|F|(42),) = Y i-y^^^^'-' (2L + {{12),\Vm),) . (13) 

L [ ji 32 L 

Conversely, we have 
((12)^1^1(34)^) = Y {-r^^'^"^ (2A + 1) ^ ^ ((I3),|y|(42),) . (14) 

A J4 is A 
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It is useful to introduce here some notations we will be using throughout the 
rest of this paper. First, we note that, by definition, the scalar product of two 
tensor operators of rank k, U^''^ and V^''\ is introduced as 



^U(k) y = ^ Jjik) yik) (^5) 

q 

In turn, the definition of the spherical tensor product of two tensor operators, 
^(^i) and B^''^\ is 



Second, we note that the angular degrees of freedom will be described in terms 
of the unnormalized spherical harmonics, defined in terms of their normalized 
counterparts such that ([23], Eq. 2.5.31) 



C^'\r) = ^Y,,{f). (17) 
The unnormalized spherical harmonics satisfy the orthogonahty conditions 



For completeness, we list here useful properties of the unnormalized spherical 
harmonics: 

i) cW(f) = 1. 

ii) Cf)*(f) = (-)«cl5(f). 

iii) {C^^\r) C'^^\p)^ = Pfe(cosu;), where u is the angle between f and p. 

iv) [C('=i)(f) ® C('=2)(f)]^''^ = (A;iOA;2 | kO)C^''^{f). 

v) C(f^)(f) C(^^^)(^) = Efc, (f) (A;iOA;2 | kO) {k,q,k2q2 \kq). 

Finally, in the interest of brevity, throughout the remainder of this paper 
we will follow the conventions, and will be referring directly to the equation 
numbers, found in Ref. 



3 Matrix elements calculation using harmonic-oscillator wave func- 
tions 



Provided that the single-particle wave functions are harmonic-oscillator wave 
functions, one can use Talmi transformation [IS] via the Moshinsky trans- 
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formation brackets [20] to exactly factorize a product of two single-particle 
wave functions into a part which depends only on the center-of-mass degrees 
of freedom, and a part which depends only on the relative-motion degrees of 
freedom. For completeness, we will review now the technical details of this 
procedure. 

3.1 Moshinsky transformation brackets 

Consider a two-particle system in a harmonic oscillator potential. We shall 
characterize the two particles by their coordinates and quantum numbers. For 
the purpose of our discussion, we introduce two system of coordinates: 

• laboratory frame, where the two particles are described by their coordinates 
with respect to the center of the potential well, fi and f2, and corresponding 
radial, ni and n2, and orbital quantum numbers, h and l2- 

• center-of-mass frame, where the system is characterized by the relative co- 
ordinate f and the coordinate R of the center of mass of the two particles, 
defined as ^ 

f = (ri - ra) , R = + fg) ; (19) 

the radial and orbital quantum numbers n, I will correspond to the relative 
motion, and A^, L to that of the center of mass. 

The eigenkets in the two coordinate systems may be written as follows: 

• laboratory frame 

\n1n2 (/i/2)A/i) {Iimil2m2 I A/i)|ni/imi) \n2l2m2) ; (20) 

mi 7712 

• center-of-mass frame 

\nN (IL) Xjj,) =Y, {ImLM \ Xfi) \nlm) \NLM) . (21) 

mM 

The Moshinsky transformation brackets are defined as the expansion coeffi- 
cients of the eigenket ( l20l) in a series of eigenkets ( 1211) . We have [20] 

\n1n2 ihk) A/i) = J2 {nlNL\\nihn2l2 A) \nN {IL) Xfi) . (22) 

nlNL 

This transformation is independent of the magnetic quantum number fi, and 
the transformation bracket vanishes for all combinations of its parameters 
which do not satisfy the total angular momentum 

A = fi + fa = i + L, (23) 
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and energy 



£;oc(2ni + Zi + 3/2) + (2^2 + /2 + 3/2) 
= {2n + l + 3/2) + {2N + L + 3/2) . 



(24) 



conservation laws. Therefore, the transformation bracket vanishes for all com- 
binations of its parameters which do not satisfy the energy condition flMj) . 
and any summations over A will be restricted by the corresponding triangle 
relations: 



\l-L\<X<l + L. 



(25) 



We calculate the two-body matrix element in a harmonic oscillator single- 
particle basis using the Moshinsky transformation brackets, and obtain 

{mn, (1,1)^^(1,1)^^; JMj I V I nsn,(kl)^^(lS,:, JMj) 

= 3ihhu AA' (^1/1^2/2 A' I n'/' A^'L' A')(n3/3n4/4 A I n/7VL A) 

nlNLX;n'l'N'L' X' 



(26) 



X {n'N' (l'L')y (I l)s' ; JMj I y I niV (lL)x (^ ^)s ; JMj) . 

For a particular choice of the potential, V, all we have to do is to calculate 
the matrix element 



1 


h 


1 
2 


Ji 




^3 \ J3 


X Y ' 


k 


1 
2 




> < 


^4 \ 34 ' 


S,S'=0 












A' 


S' 


J 




X S J 



{n'N' (l'L')y (1 \)s> ■,JMj\V\nN (IL), (i 1)5 ; JMj) . 



(27) 



3.2 Central interaction. 



For the case of the central interaction, the potential depends only on the 
magnitude of the vector f, and not on its angular degrees of freedom. We 
apply the Wigner-Eckart theorem for the case of a zero-rank tensor (|23j. 
Eq. 5.4.1a) 



{n'N' (l'L')^,(\ \)^ - JMj \V,\nN (IL)^(\ i)^; JMj) 



(21 



-J {n'N' (l'L')^,(\ \)^ - JMj II K II nN (IL)^(\ i)^; JMj) 
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where we have introduced the notation j = \J1J + 1. The reduced matrix 
element in Eq. (128|) is evaluated using ([23], Eq. 7.1.7). We have 

(nW {l'L')^,{\ 1)^,; JMj \\ V, \\ nN {IL)^{\ \)^- JMj) (29) 

6ss' (-)^'+^+^+o (2J + 1) ,) ^' ^ ^ ^n'N' {I'L'),, \\ V, \\ nN (IL),) 

J A 



We use ([23], Eq. 7.1.7) one more time for the reduced matrix element 



n'N' {I'L')^, II Vc II nN (IL)^) (30) 

W5ll'(-)''+'^+'+°AA' <( ^ ^ ^ }{n'l' II K \\nl). 

A / 



The Wigner 6j symbols in Eqs. ( 12^150]) are calculated using ([23], Eq. 6.3.2). 
From the Wigner-Eckart theorem ([23], Eq. 5.4.1a) we have: 



{n' I'mi \Vc\n Imi) = 5w j (n' /' || K || n I) (31) 
= 6n> nM[V,]{n'l-nl), 

where we have introduced the notation 

nM[V]{n'l' ■,nl)= p dr] HOn'i'{r) V{r) HOni{r) . (32) 
Jo 

Therefore, Eq. 0281) becomes 



{n'N' {l'L')^,{\ \)^ - JMj \V,\nN {IL)^{\ \)^- JMj) (33) 
= ^ss'^NN'^LL'h'^xx' TZM\Vc]{n'l ; nl) . 



3.3 Spin-spin interaction 



It is convenient to carry out the calculation of the spin-spin interaction matrix 
element by using the m-representation approach outlined above. We have: 



9 



{n'N' (l'L')^,[\ \)^ - JMj \V,a^■a2\ nN (IL)^{\ \) ^- JMj) 
J2 {I' mv L' M' I A' fi') (A' fi' S' nis' \ J Mj) 

mil M' rrigi ^' 

X {ImiLM \ \^){\^iSms \ J Mj) 

X {n'N' I'mr L'M' S'ms> | K ai ■ as | nN nN Imi LM Smg) , 
where the matrix element in the last equation is given by 



(34) 



{n'N' I'mv L'M' S'nis' | K f^i ■ | nN nN Inii LM Sms) 

= ^NN'hL'hlM'^W^mimi, {n'l \ K | nl) {{\\)s' I CTl ■ 0-2 I {\\) s) ■ 

The spin-dependent factor is obtained as 



(35) 



{(k\)s' I CTi ■ 0-2 I {\\)s) 



S{S+l)-- 



^SS' 



(36) 



where we used the fact that Sj = ^ai. Therefore, we obtain 



{n'N' {l'L')^,{\ \)^ - JMj I ai ■ (T2 I nN {IL)^{\ \) ^- JMj) (37) 



Snn'Sll'Smm'Su'Sss'^xx' 2 



S{S + 1) 



nM[Vs]{n'l ; nl) . 



3.4 Tensor interaction 



We begin the calculation of the tensor interaction matrix element by using 
(p], Eq. 7.1.6) 



{n'N' {l'L')^,{\ \)^ - JMj I VtS^2 \ nN {IL)^{\ \)^- JMj) 
' J S' X' 



(38) 



^+^'+^ ^ ^ E i^'N' {l'L')y II V, C(^) II n"N" (IL),) 

2 X S \ n"N" 



y,{n"N"a\)s' ||[ai®a2](2) || nN [\\)s) , 



where 



^12 = 3 (ai ■ f) (a2 ■ f) - ^1-^2 = 76 (C(2)(f)0[ai®fr2](2n . (39) 
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The reduced matrix elements in Eq. (l38l) are evaluated as follows: 
(i) We use (|23], Eq. 7.1.5) - for S{S') = 0, 1, to obtain 



1 i 1 

2 2 

1 i 1 

2 2 

S' S2 



\a II ^1 II 



with (i II ai = ^/6 (see Ref. [23], Eq. 5.4.4). 



(ii) Next, we use ([23], Eq. 7.1.7) to obtain 



{n'N' {l'L')y II Vt II niV (/L),) 

/' A' L 



^NN'^LL' 



\l'+L+X+2 



A / 2 



(n' I' II II n I) 



(iu) We also have 



{n' I' II II n /) = (/' II II /) 'JZM[Vt]{n'l' ; n/) 



where (/' || C^'^) || I) is given by Eq. flRT]) 



We obtain 



(n'iV' (/'L');,,(i i)^,; JMj | | riN {IL)^{\ \) ^- JMj) 

= Snn'Sll'Sss'Ssi 2^30 (-)''+^+-^+^ XX'l (/0 20 | /'O) 




X <^ }{ }nM[Vt]{n'l' ; nl) . 
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3.5 Spin-orbit interaction 



Similarly to the case of the tensor interaction, we first use ([23], Eq. 7.1.6) 
{n'N' {l'L')^,{\ \)^- JMj I Vis (^ S) I nN {IL)^{\ \)^- JMj) (44) 



\\+S'+J 



J S' X' 
1 A S 



X Yl (^'^' i^'L')y II ViJ II n"N" {lL)y){n"N" {\\)s' ||S || nN {\\)s) 

n"N" 



with the reduced matrix elements calculated as: 
(i) We use (|23], Eq. 5.4.3) - for ^ = 0, 1, to obtain 



{i-2i)s' lis II C2i)s) = Sss'^SiS + 1)(25 + 1) = 76 6siSss' • (45) 
(ii) Next, we use ([22], Eq. 7.1.7) to obtain 



{n'N'{l'L')y II Vis£\\ nN{lL)x) (46) 

6nn'6ll' {-y'+^+'+' XX' 1^^' ^'^ ^ i^n'l' \\ ViJ \\ nl) 

X I 1 



(iii) Finally, we use (/ || i \\ V) = 5u> I + 1) ([23], Eq. 5.4.3) to obtain 

{n' l'\\ViJ\\nl) = 6u' 1 ^Jl{l + 1) UM [Vi^] {n'l ; nl) . (47) 
Collecting terms, we find 

{n'N' {l'L')^,{\ \)^ - JMj \Vis{£-S)\ nN {IL)^{\ |)^; JMj) (48) 
= 5NN'hL'5u'5ss'5si (-)'+^+-^ m v//(/ + l) 

X 

3.6 i"^ Interaction. 



The matrix element of the i"^ interaction can be easily calculated when the 
matrix element is written in the m representation, i.e. 



Ijsx'U 


\ I X' l] 








^7ZM[Vis]{n'l 


\i xsj] 


[X I ij 
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{n'N' I'mv L'M' S'ms' \Vi2f\ nN Imi LM Sms) (49) 

= ^NN' hv hiM' 6m,mi, ^SS' ^msnis, + 1) TlM[Vi2]in'l ; Til) . 

We use the transformation 

{n'N' {l'L')^,{\ \)^ - JMj I Vi2f I nN {IL)^{\ |)^; JMj) (50) 
^ (/' mi, L' M' I A' /i') (A' fi' S' nig, \ J Mj) 

mi/M'mg/fj.' 

X J2 {Irrii LM \ \fi){\iJ,Sms \ J Mj) 

niiAImsfj, 

X {n'N' I' nil, L'M' S'nig, \V12f \ nN Imi LM Sms) , 

substitute Eq. (149|) . and employ the orthogonahty properties of the Clebsch- 
Gordon coefficients, to obtain 



{n'N' {l'L')^,{\ \)^ - JMj I Vi2e \ nN {IL)^{\ 1)^; JMj) 
= Sxx'6NN'hL'Sii'6ss' l{l + 1) nM[Vi2]{n'l ; nl) . 



(51) 



3. 7 P' {o\ ■ 0-2) interaction 



The derivation of the matrix element for the P' {a\ ■ (T2) interaction follows 
closely the calculation of the matrix element corresponding to the i'^ inter- 
action. We first perform a transformation to the m representation where the 
calculation of the matrix element is particularly simple, i.e. 



{n'N' {l'L')^,C2 Dsn I Vi2sf (^i ■ ^2) | nN {IL)^{\ \) ^- JMj) (52) 

{X'l^' S'ms, I JMj){X^Sms \ J Mj) 

X {n'N' \' jj! S'ms' \ Vi2s ^ (cri ■ cr2^ \ nN Xjj, Sms) , 



which gives: 



{n'N' {l'L')^,{\ \)^- JMj I Vi2se [a, ■ as) | nN {IL)^{\ \)^- JMj) (53) 

S{S + l)-\ 



<^AA' ^SS' ^NN' ^LL' 2 



/(/ + 1) nM[Vi2s\{n'l; nl) 
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3.8 Quadrupole spin-orbit interaction 



For the purpose of generalizing this matrix-element calculation to the case of 
an arbitrary set of single-particle wave functions, it is convenient to approach 
this calculation by first noting that 



j=0 



l(i) 



(0) 



(54) 



A close inspection of the above equation shows that some of the pieces of the 
quadrupole spin-orbit interaction may be incorporated into the calculation of 
the other previous interactions involving the relative orbital angular momen- 
tum operator, i. 



(i) Case j = 0. 



[i ® ef^ ® [s s] 



(0) ] (0) ^ 1 
6 



f {3 + a^■ , 



(55) 



where we have used {af^ = afy = af^ = 1) 



(56) 



Therefore, we can introduce the modified radial amplitudes of the P and 
i'^ ^(Ti ■ (72) interactions: 



V12— V12 + - Vis2 , 



Vl2ss — Vl2s 



V, 



ls2 ■ 



(57) 
(58) 



(ii) Case j = 1. 



(1) 



S® S 



(1) 1(0) 



■ S 



(59) 



where we have used the definition of the angular momentum quantum opera- 
tor, J X J = i J. Similarly, we introduce the modified radial amplitude for the 
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spin-orbital interaction, as: 



Vis — Vis — 2 



(60) 



(iii) Case j = 2. 



(0) 



75 



[i(S)if^ ^[ai^a2f'>f'^ . (61) 



We notice that the only component of the (j! ■ interaction that we have 
not addressed yet is the one corresponding to j = 2. 

To calculate the corresponding matrix element we first use ([23], Eq. 7.1.6) 

{n'N' {l'L')^,{\ \)^- JMjIViJ [[i ® i]^'^ ® [S ® S]^^')] \nN (/L),(ll)^; JMj) 
J S' A' 



^_^x+s'+j+,^ " - K Y: {n'N'{l'L')y II Vi,2 [i^if^ II n"N"{lL),) 

j \ S \ n"N" 



X {n"N"C^l)s' \\[S®Sf^ llniV(ii)^), 
with the reduced matrix elements calculated as: 



(62) 



(i) We use ([23], Eq. 7.1.1), to obtain 



((li)^' ||[S®S](^') II (ll)s) = 6ss'Ssi6{-yj 



1 1 J 
1 1 1 



(63) 



(64) 



(ii) Next, we use ([23], Eq. 7.1.7) to obtain 

{n'N' {l'L')y II Vi,2 II nN (IL)^) 

= Snn'Sll' (-)''+^+^+^- AA' r' ^' ^ ^ {n'l'\\Vis2 [i ® i]'''^ \\nl) 

[a / J 

with ([23], Eq. 7.1.1) 

{{n'l'\\Vi,2[i^if^ II nl) 

{-y 3 l{l + 1)(2/ + 1)1 l^^^ \ nM[Vis2]{n'l ; nl) . (65) 

[ill] 
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Collecting terms, we find 



2 2' S' 



JMj\Vls2 



n'N'{l'L'),,{ 
= Snn'^ll'Su'^ss'^si 6(— ) 



l+L+J+l 



/(/ + 1)(2/ + 1) 



(66) 




nM[Vis2]in'l ; nl) . 



4 Two-body matrix elements calculation 



Unlike the case of harmonic-oscillator single-particle basis set, the calcula- 
tion (and storage) of the two-body matrix elements using single-particle wave 
functions that can be represented as linear combinations of harmonic-oscillator 
wave functions is performed more efficiently in ph coupling. The j»/i-coupled 
matrix-element calculation is based on the following two lemmas: 

Lemma 1: Provided that the two-body potential can be factorized into parts 
depending only on the r*i or r*2 coordinates, respectively, 

V{f,,f2) = {u^'\n)QV^'\f2)) , (67) 

then ph- coupled matrix elements of the two-body interaction are given by 

((13),| {U^'\n) V^'\f2)) 1(42),) (68) 
-(1 II II 3) (2 II II 4) 4a. 



2A + 

Using the definition of the scalar product of two tensor operators of rank k, 
Eq. f fTSj) . we calculate the matrix elements of the U^'^^ {fi)vl''J {f2) operator in 
the m-representation, using the Wigner-Eckart theorem 

(12|[/f (fi)vd?(r-'2)|34) = ^-1- (jimua - m, \ kq){l \\ U^^^ \\ 3) 

k 

{32 m2 j4 -m^lk - q){2 \\ V^^^ \\ 4) , 

(69) 

where we have also introduced the notation k = \/2k + 1. Then, Eq. ( 169|) is 
obtained using the orthonormality of the Clebsch-Gordon coefficients, together 
with the definition of ph-coup\ed matrix elements, Eq. (ITTi) . 



C_')i4-m4 

X 
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We note the identity: 



(0) 



(70) 



where the spherical tensor product of two tensor operators, A^'^^^ and B^''^\ 
was defined in Eq. (fT6|) . 

Lemma 2: Provided that the spatial part of the two-body interaction has the 
form V{r) C^''\f), with k a positive integer or zero, then the variables ri and 
r*2 can be separated in the sense that [24^ 

r-V{r)C^'\f) i^2-^i-^ (^^i + y^^^ + ^) (fc,0fc2 | kO) (71) 



ki,k2 



(fc) 



where a is a positive integer exponent, and r = fi — r2 indicates the relative- 
motion coordinate. We have also introduced the notations 



2 r°°, 



TX Jo 



[P^ dp] f a K (p) Jfci {pri ) jk^ {pr2 ) 



(72) 



and 



VaKip) = / [r"^ dr]r°'V{r)j^{pr) . 
Jo 



(73) 



To prove lemma 2, we first introduce the asymmetric Fourier transform of the 
operator V{r) C^^\f) 



v(»''\p) = —1- /dVr"l^(r) C^''\r)e'P- 



(74) 



and, conversely. 



r»V{r) C'^''\f) = d?p i^^\p)e~'^'' . 



(75) 



We use the plane wave expansion in terms of unnormalized spherical harmon- 
ics, 



$:i'(2/ + l)jHpr) (C«(f)0C«(p)) , 



(76) 



together with the orthogonality conditions (lT8l) . and carry out the angular 
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part of the integral fl74p . Then, v'^^\'p) becomes 

v^'^'Kp) C^'\P) v.M ■ (77) 

From Eq. ( 1751) . we can write 

r-V{r) CW(f) = ^/d^P v^M C^^^p) e^^" ^ (78) 

The angular part of the last integral can be carried out explicitly using the 
definition of the relative coordinate, f = ri — r2i and by applying Eq. fl76l) 
twice for exp(ip ■ ri) and exp(— ip ■ r2)i respectively. The following integral is 
evaluated using the above properties of the unnormalized spherical harmonics: 

j dfi, {Cf^^\p)Y [C^i'\p)r C^^\V) = {hOhO I A;0) {hq, hq2 \ kq) . 

(79) 

To finalize our proof, we use the definition of the spherical tensor product of 
rank k obtained from the tensor operators C^'^'^\fi) and C'^^'^\f2). 



5 Matrix elements calculation using arbitrary single-particle wave 
functions 

It is convenient to evaluate two-body matrix elements in ph angular momen- 
tum coupling using lemma 1. Then, according to lemma 2, ph matrix elements 
factorize in two parts, which depend on the coordinates of the first and the 
second particle, respectively. Hence, the appropriate angular-momentum cou- 
pling for each single-particle wave function |a) is provided in the {Isj) coupling, 
with the individual orbital angular momentum and spin Sq are coupled to a 
total angular momentum ja- For an operator O'^^'^ which depends only on the 
orbital angular momentum components, in the (/sj)-coupling scheme we have: 

i) (see [23], 7.1.7) 

II O^^'^ II {h\),,) = F{hn; I2J2; k) {h II O'-''' II I2) , (80) 

with 



F{liji]l232\k) 



-TT+fc 
2 J1J2 



h jl \ 

32 h k 



(81) 
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ii) (see [23], 7.1.5 and 5.4.4) 



{{hi),, II ® a ^"'^ II (/4),,) = G{hn-h32\hh) {h II O^'^^ II k) 



(82) 



with 



G{hjl]l2j2]hk2) = VQjlj2k2 



h h h 

1 i 1 

2 2 

ji J2 A;2 



(83) 



We can now proceed to discussing the calculation of the two-body matrix ele- 
ments of the Argonne vis potential: We will first evaluate the radial part of the 
matrix element, which can be performed independently of the way we handle 
the angular and spin/isospin degrees of freedom of the interaction. Next, we 
will discuss the angular and spin part of the 7 basic operators corresponding 
to the Argonne vis potential. 



5.1 Radial part of the two-body matrix element 



The radial part of the two-body matrix elements is defined as 

2 roo 



TT Jo 



[p dp\vcK{p) 



X 



roo 

/ [r^drl]7^„,^,(rl 

/ [rl dra] 7^ 
Jo 



Jk2{pr2) 



B4) 



where VaK,{p) was defined in Eq. (1731) and the operators (9^^'^^ are operators 
acting only on the radial part of the single-particle wave functions. The func- 
tional form of the radial wave functions in Eq. flH^ is not restricted to that of 
a harmonic-oscillator wave function. Therefore, at this point we need to care- 
fully consider what is the most efficient numerical strategy to computing fl8^ 
for the particular radial form of the wave functions under consideration. For 
illustrative (and testing) purposes, in the following we will confine our discus- 
sion to the case of a basesis of single-particle wave functions defined as linear 
combinations of harmonic-oscillator wave functions [see Eq. ([1])], in which case 
some of the integrals in (18^ can be efficiently performed using Gauss- Hermite 
quadrature formulas. 
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Treating separately the integral 



^5) 



we first change variables such that x'^ = Ti {y/2/h) with q' = p {b/y/2), and we 
recall that (see discussion surrounding Eq. (E]) 



Next, we perform the harmonic-oscillator expansion 



(86) 



-^"I'l (75 1 



(87) 



Therefore, the integral (185|) becomes 



EAnki 



(88) 



where we have introduced the notation, T-tOneiq) = {—)"'HOne{q), which de- 
notes the Fourier transform of the harmonic-oscillator wave function, TiOne, 
and the coefficients A'^^l^^^i^ are calculated as 



f°° 2 
•J 



nOnk,{x'i). (89) 



This integral is conveniently performed using Gaussian quadratures. Similarly, 
we have: 



Jo 



/\mk2 



-no.^,,{q') 



(90) 



where the coefficients A^^j^^n^u given as 



^00 r "I 

KlUu = / [Xldx'ji?„,j^x'2) 0(^)/?n,u(^x'2) ^0™.,(X'2). (91) 



Using Equations (188|) and (pOl) . the radial part of the two-body matrix elements 
becomes 
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K 



,aK;kik2 

nih ri2h ;"3'3 '^4U ' 



EAuki Amk2 



[q'^ dq']v^^{^c^)HOnuA^)HOmk2{q') ■ (92) 



5.2 Central interaction 



The central interaction depends only on the magnitude of the relative dis- 
tance r between the two particles, i.e. 



Vc = V,{r)C^'^\f). (93) 
The strategy to calculating the ph matrix element 

is emblematic for the rest of this section: (For simplicity, we will not concerns 
ourselves here with the radial degrees of freedom. The radial part of the matrix 
elements is calculated using the approach discussed in the previous section.) 
First, we use lemma 2 and separate Vc into tensor-operator components which 
depend on either ri or r2. We have 

Vc = 5](2A: + l)M('='='°°)(ri,r2) {C^^\r\) Q C^''\f^)) . (94) 

Second, we use lemma 1 to calculate the ph matrix element of the correspond- 
ing interaction. For the central interaction, this procedure leads to 

X II C^'\n) II (/3|),3) {{hlh II C^'\r2) II mW > (95) 

which gives 

{ihlUhlh, ; A/i I Vc I {l4hM)h ; V) = i-y^^^^^' n(^^^°°)(ri,r2) 
X F{hn; /3J3; A) (/i II C(^) II h) F{l2j2; Uja; A) {k \\ C^") || U) . (96) 

5.3 Spin-spin interaction 

By definition, the spin-spin interaction is introduced as 

Vs = Vs{r)a,-a2, (97) 
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where we can write 

We perform the recouphng 



0-2 



{91 







(0)" 









(0) 



(99) 



cr2 



(0 



and obtain 



((^i|)ii(^3|)i3 ; A/i I I (^4|)i4(^2|)j2 ; V) 



(100) 



2A + 



X G(/iji; /3J3; A;A) (/i || C^^') || /g) ^(/sjs; ^4j4; kX) {h \\ C^^) || U) , 
with |A - 1| < A; < A + 1 . 



5.y^ Tensor interaction 

The tensor interaction is defined as: 

Vt = Vt{r) Su , 
where the operator 5*12 defined in Eq. (139|) is written now as 



(101) 



s 



12 



(0) 



Using a recouphng scheme similar to Eq. f l99|) . 



1 1 k 



(Ti 0-2 



(2) 



(0) 



0-1 



(fc) 



(102) 







C('^)(r2)0 



0-2 



(103) 
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we obtain 



X 76 ^ {2h + 1) Gikn; /3J3; h\) {h II II /g) E i"'"""' (2^2 + 1) 

fcl fc2 



(104) 



X ^('^^'^^^^o^^^^^^^^^^Q^^Q I 2 0) <( ^ ^ Gikj2;kj4;k2X){l2 II || /4) • 

1 1 A 



5.5 Spin-orbit interaction 



The spin-orbit interaction is given by: 



(105) 



where the orbital angular momentum operator and the total spin operator are 
defined as 



i f X (pi - P2) , and S = ^ {<^i + ^2) 



(106) 



respectively. Using the tensor operator properties: 



C^\f) , and (f X yU = -i V2 [x ^ y^^ , (107) 



we can write 



£= — = r 
V2 



C(i)(f)®(v2- Vi) 



(1) 



(108) 



Then, we use lemma 2, and obtain: 



fclfc2 

(2A:i + l)(2A;2 + l) 



X 



3^2 



(fciOfczO I 10)M('=^^^2'")(ri,r2) 



C^'^l^fi) ®C(*^2)(f2) 



(109) 

1(1) 
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Recoupling the C^^^)(fi) and C^'''^\r2) operators with the appropriate gradient 
operators, we obtain: 



Visir) i = i"^-"'-' i2h + 1) i2k2 + 1) {h0k20\l 0) 

fclfc2 



k \ ki k2 1 



? ^ 1 1 1 k^ 



(fcifca.ll) 



(ri,r2) 



(110) 



+ 



1(1)' 



Next, using the triangle conditions for the angular momentum arguments of 
the Wigner 6j symbol, we derive the conditions: \ki — k2\ < 1 < ki + k2, 
\k2 — l\<k<k2 + l and |/ci — 1| < k < ki + 1. These conditions together with 
the condition that ki + /c2 is odd, derived from the Clebsch-Gordan coefficient 
in ( IllOp . can be satisfied only if k equals either ki or k2- Therefore, we can 
write Eq. ffTTOj) as 



Visir) 



— y 



•fcl+fc2 + l 



X 



(2A;i + l)(2A;2 + l) (/c2 1 | k^O) 



ki k2 1 
1 1 ki 



(fcl 



(1) 



1(1) 



+ U 



(fcafci,!!) 



(^■1,^2) 



(1)' 



(111) 



with ki + k2 + 1 even. 



Taking the inner products of Eq. f illip with S, and using the following iden- 
tities valid for operators A^*^^ and B^^^ that commute with all components cxj 



^(1) ^{2) 



^1 =(-) -T 



5, 



(^) 

(2) 



(1) 



fe+1 

k 



a2]={-r'^{A\'^^Q 



^{2) f72 



(112) 
(113) 



the operators in Eq. fillip give the following contributions to the ph-conp\ed 
matrix element ((/i|)ji(^3|)i3 ; A/x | Vls \ {kl)jiil2l)j2 V): 
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(1) 



S 



(114) 



JkiX 



.)A+j2+J4 + l ^3 



A(2A + 1) 2 



X 



G(/ui;/3j3;AA) {h II C^'^ II k) F{l2j2;kjA;X) CG2iXh,k2X) 
-F(/iji;/3j3;A) (/i II II /3) G{kj2;hK,\\) CG2{\k2,k2\) 



and 



(_)fcl+i2+i4 + l ^ 



(1) 



S 



(115) 



A(2A+1) 2 



X 



Kx G{kn-hj3;ki\) GGi{k,X,k,ki) F(/2j2; /4J4; A) (^2 || || Z^) 
-4,A F{hj,;hj3;\) CG^i\k2,\\) G{l2j2;kj4;k2\) {h \\ G^'^^ \\ k) 



and 



(_)fc2+i2+i4 + l y/3 

A(2A + 1) ^ 



(1) 



S 



(116) 



X 



Gihji;kj3;k2X) {k II C^'^^) II h) F{1232-M3aA) CG2(A;2A,AA) 
-42A F{hji;hj3;X) {h II C(^) II /3) G{l232;h34;kiX) GG2{\ki,k^k^) 



and 



C(^'^Hn)0ViJ^''^®C'(^'H^2) 



(1) 



5 



(117) 



5u- 



kiX 



A(2A + 1) 2 



X 



G{hjuhj3;XX) GG,{k2X,k2X) F(/2j2; /4J4; A) (^2 || || z^) 
-i^(/iJi;/3j3;A) CGi(A;2A,A;2A) G(/2j2; /4J4; AA) {k \\ || h) 
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where we have introduced the notations 



CG^{hk2jil2) = {h\\u^''''''^''\n,r2) C('^)(n)®V(i) ^''^ 11/3), (118) 



'(2) 



('2) 



h) . (119) 



5.6 i'^ interaction 



The i"^ interaction is defined as 



Vl2 = Vi2{r)f 



(120) 



where we can write 



f = -y/3 [e®i] 



(0) 



(121) 



It is useful to discuss this as a particular case of the more general operator 



(i) 



>® (pi -P2)]^^^ ® [r ® (pi -P2)]^^^ 



(122) 



We begin by using the definition of the spherical components of an arbitrary 
vector, a, i.e. ([23], 5.9.4) 



1 / 

a±i = =F--/| (ctx ± 1 a^) ; oq = 



(123) 



Then, we can show that 



>„,,p„] = ih (-)™ 



(124) 



and 



[ (ri - rs)^, (pi - P2)n] =2ih{-r6n,. 



(125) 
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We can change the couphng scheme and combine the two rs, in Eq. f ll22p . 
into a single tensor operator dependent on the relative-coordinate unit vector: 



+ 6 i r 



(1) 

1 1 1 
1 1 1 

K k' j 



(1) 



r®{pi-p2)r' (126) 

J m 



(i) 



1 1 1 



1 1 J 



Because of the cross product property, Pm >^ Pn = 0, 1^ and n' cannot be equal 
to 1. (The parameter k is also restricted to even values only because of the 
Clebsch-Gordan coefficient in fll26p .) Then, we can write 



0') 3 



2 ^ 



X 



3 r 



1 1 1 
1 1 1 

K k' j 

C('^)(f) ® 

1 1 1 
1 1 J 



> (10 10 I kO) 

'Vi - V2) ® (Vi - V2 

ca)(f)®(vi- V2)'^'^ 



(127) 



In the particular case of the £^ interaction, the rank j is equal to 0, and we 
can use the symmetry properties of the Wigner 6j and 9j symbols (|23], 6.4.14 
and 6.3.2) to obtain the tensor-product form of £^ as 



2 K=0,2 \ 1 1 K 



;io 10 1 ko) 



r) ® Vi - V2 



(0) 



IVi- V2)®(Vi- V2j 



(0) 



:i28) 
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We expand the unnormalized spherical harmonics C^'^\r) and C^^\r) using 
lemma 2, to obtain 



1 1 1 1 



X 



K=0,2 ^ I 1 K 1 



(129) 



(1010 I kO) Ei'''"'''~''(2^i + l)(2A;2 + l)(A;iOA;2 | kO) 

fcifc2 



X M('=i'=2'"2)(ri,r2) 



r2 



(Vi- V2)®(Vi- V2' 



(«)' 



(0) 



E i''"''^~'(2A;i + l)(2A;2 + l)(A;iOA;2 | 10) 

v3 



X 



(0) 



(1) 



(0)- 



The operators in Eq. (11291) give rise to contributions to the matrix element 
((^i|)ii(^3|)j3 ! -^/^ I ^L2 I {h'^)j4{h^)j2 ; -^At) similar to the contributions of the 
i'^ ((Ti ■ 0-2) and quadrupole spin-orbit operators discussed in the subsequent 
sections. Therefore, it is useful to base the calculation of the P matrix element 
on several general results. 

First, let us discuss the case of the matrix elements involving only one gradient 
operator in Eq. (I129p . which are very similar to those discussed in the case of 
the spin-orbit interaction. We have the general identities: 



h k2 1 



(i) 



(130) 



^V3 E ^ 



1 J k 



C^'''\ri)®\C^^'\r2)®V2 



(fc) 



and 



C('=i)(fi) ® C(^'2)(f2 



(1) 



® Vi 

k2 ki 1 
1 J k 



:i3ii 



ri) O Vi 



(fc) 



r2 



U) 



For the i interaction we are only interested in j = 0. We obtain: 



(_)A+j2+i4 + l 

^ ^fc i A — 

' A(2A + 1) 



(1) 



® V2 



(0) 



(132) 



F{hji;hj3;X){li II C^"'^ II k)F{l2j2;lAj4;X)CG2{Xk2,k2X) 
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and 



(1) 



Vi 



(0) 



(133) 



S. 



k2\- 



A(2A + 1) 



-F(/ui;/3j3;A)CGi(A;iA,fciA)F(/2j2;/4j4;A)(/2 || C^^) || ^4) • 



Next, we multiply the operator expressions in the first sum in Eq. fll29p to 
obtain 



{K) 



® 











(0) 


) (Vi - 


- V2j_ 








(«) 


(0)" 




V2 ® V2 






in) 












Vi 


® V2 





(134) 











Vi ® Vi 





+ 



Now, we have to change the coupling and separate the operators depending 
on the coordinates of the first particle, from the operators depending on the 
coordinates of the second particle. Hence, the matrix elements involving two 
gradient operators in Eq. fll29p . are calculated using the following identities: 



V2 5 


?> V2 




K 

4 


C('=^)(fi) 









(i) 



(135) 



V2®V2 



(«') 



iJ) 



and 



Vi ® V2 



ki k2 K, 
1 1 k' 
k k' j 



iJ) 



(136) 



(fe') 



(i) 



and 







(«')" 




Vi ® Vi 





(i) 



j k 



Vi®Vi 



(k) 



(137) 

liJ) 
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Again, for the interaction, we are only interested in the case j = 0. We 
obtain: 











V2 ® V2 





(0) 



^^^>^ X{2\ + 1) -^(^i^'i'^3j3;A)(/i II II /a) 
X F{k32] hji, X)CGG2{Xk2K, /c2kA) 



(138) 



and 



2A + 1 1 1 A 











Vl ® V2 





(0) 



(139) 



X F(/iji; /3J3; A)F(/2j2; /4J4; A)CGCG'(A;iA;2K, k,Xk2X) 



and 



'C(*^i)(fi) ®C(*^2)(f2) 
('_^A+i2+i4+i 











Vl ® Vl 





(0) 



(140) 



— F(/iji; /3J3; A)CGG'i(A;iAk, /cikA) 

A (zA + Ij 

X F(/2J2;/4J4;A)(/2||CW \\ u) . 



Here we have introduced the notations 



CGG,{k,k2n,klX) = {h II u^''''''^^\n,r2) \G^'\h) ® [V(i) ® V{i)] 



1(0 



GGG2{kik2K, klX) = {I2 II u('=^'=^'2«)(ri, r2) fc('=)(^2) ® [V{2) ® V(2)] 



1(0 



(A) 



(A) 



II h) 

(141) 

II ^4) 
(142) 



and 



GGCG{k^k2K, /1/2A1A2) = {h II n(^^^'2'^)(ri) C^'^^^^i) ® Vi ^''^ II ^3) (143) 



X (/2 II ^^'^''"^^2) [C(^^^(r2)® Vi 

5. 7 £^ (^cTi ■ (72^ interaction 

The £^ ^ai ■ (T2) interaction is given by: 



{A2) 



h). 
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L2S 



(0) 



(144) 



Since the £^ (ai ■ asj interaction differs from the i"^ interaction only through 
the spin part, (ai ■ as), we can use Eq. ( 1129^ and add the corresponding spin 
interaction. We have: 



ft=0,2 I 1 K 1 



(145) 



X 



(Vi- V2)® (Vi- V2) 



(0) 



(CTI ■ 0-2) 



V3 



fclfc2 



^2— fcl — 1 



(2A:i + 1)(2A;2 + 1) (A:i A;2 | 1 0) u^'''''''^^\n,r2) 



X 



C('=i)(fi) ®C('=2)(f2) 



(1) 



V2 



(0) 



(cTi ■ (T2) 



C('^l)(fi) ®C(^2)(f2) 



(1) 



1(0) 



Vi 



((Ti ■ 0-2) 



Accordingly, using Eqs. fll38H140"ll and Eqs. (113211133 1 ) we obtain the contri- 
butions to the jj/i matrix element {{hl)jiihl)js', ^1^ I Vl2S I {hl)j4{l2^)j2i ^fA- 



and 



|'_U+j2+j4 



[V2 ® V2 



(0) 



((Ti ■ a2) 



A;i (2A + 1) 



Gihn;hjs;k^\) ih II C^"'^ \\ h) 



X G{l2j2'-,hji\ki\) CGG2{kik2K, k2K.ki) 



(146) 



(_)A+fc2+j2+i4+l 



Vi ® V2 



1(0) 



(ai ■ 0-2) 



2A + 1 



ki k2 hi 
1 1 k 



(147) 



X ^(/iji; /3J3; A;A) G(/2j2; /4J4; ^A) CGCG{kik2K, kikk2k) 
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and 



and 



and 



Vi ® Vi 



l(«) 



1(0) 



(Tl ■ (72 



A:2 (2A 



- G{liji;kj3;k2X) CGGi{kik2K, kiK,k2) 
X G(/2j2; /4J4; A;2A) (/2 || C^"'^ \\ k) , 



(148) 



(_)A+j2+i4 



(1) 



V2 



(0) 



(7i ■ 0-2 



A;i (2A + 1) 



G(/iji;/3j3;A;iA) {h \\ C^^'^ \\ h) 
X G{l2j2;kj4;ki\) GG2{kik2,k2ki] 



(149) 



G^^'\h) ® G'-^^\f2 

(_)A+i2+i4 



(1) 



Vl 



(0) 



(Tl ■ 0-2 



k2 (2A 



- (^(/iji; /3J3; k2\) GGi{kik2, kik2) 
X G{l2j2;hu;k2\) {kWG^'^'^ \\h). 



(150) 



5.5 Quadrupole spin-orbit interaction 



The quadrupole spin-orbit interaction, Vls2, is given as the radial factor, 
Vis2{'r), multiplying the operator {i ■ S)^ defined in Eq. fl5^ . As discussed in 
Sec. [31 the j = 0, 1 components in Eq. can be incorporated by introducing 
modified radial amplitudes of the spin-orbit, £^ and £^ (^ai ■ (T2) interactions, 
see Eqs. fl57ll58|60p . The only component of the interaction that we have not 
addressed yet is the one in Eq. (ETil . corresponding to j = 2. 
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In general, we can use Eq. f ll27p and write 



X 



X 



rere'=0,2 
•k2—k\—K 



1 1 1 
1 1 1 

K k' j 



(151) 



(10 10 I kO) ^i'=2-'=^-''(2A;i + l)(2A;2 + l)(A;iOA;2 | K0)u^''''''^^^\ri,r2) 

kik2 



.vi- V2)®(vi-V2; 



1 1 1 



i'^^-'^i-i {2ki + 1)(2A:2 + 1) (A;i A;2 I 1 0) 



1 1 j I fcifc2 

X M('=^^2'i^)(ri,r2)[[c('=i)(fi) ®C('=^)(f2) 



(1) 



.Vi- V2, 



(i) 



Using Eqs. f ll35H137]) and Eqs. fll30|13ip . with j = 2, we couple C'^^\fi) with 
Vi, and C^^\f2) with V2 dependent operators, respectively. The contributions 
to the p/i-coupled matrix element {{h\)ji{h\)j3', A/x | Vls2 \ {h\)jS^2\)~j2'i ^f^) 
are: 



(_)fc+i2+i4 I fci A; 2 
V5 



(2) 



[(Ti ® a2 



I (2) 



liO) 



(152) 



2A + 1 



1 1 A 



G(/iji; /3J3; h\) {k II C^'^^) II /3) G(/2j2; Uu; kX) CGG2{kik2K, k2n'k) 



and 



C^'''\h) ® Vil ® [C('=2)(^2) ® V2I ® [f^i ® ^2]^'^ 



(2) 



(0) 



(153) 



f+j2+j4 + l \ k k' 2 

2A + 1 ^ I 1 1 A 



X ^(/iji; /3J3; A;A) G(/2j2; /4J4; k'X) CGCG{kik2^, k,kk2k') 
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and 



Vi ® Vi 



(2) 



^^2 



0"! (K) (T2 



I (2) 



1(0) 



(154) 



(_)fc2+i2+j4+i A; A;2 2 

G(iui;(3j3;*:A)CGGi(*-i*-2K,fciK'*-)G(i2j2;i4j4;*:2A) ih II G<fe> II !4) 



and 



_)fc+i2+j4 \ ki k 2 
^aT^^ 1 lA 



(2) 



® CTi (g) (72 



(2) 



(0) 



(155) 



X 



G{kju /3J3; A;iA) (/i II C^'^^ II /3) G(/2j2; /4J4; kX) CG2{hk2, k2k) 



and 



C('^^)(fi)® vi ®c('='H^2) 



fc2+j2+j4 _ I A; A;2 2 



(2) 



O"! ® (72 



(2) 



(0) 



(156) 



^5 

2A + 1 ^11;, 

X ^(/iji; /3J3; A;A) CG.ik^k^, k,k) G{kj2; hu; h\) {k \\ C^''^ \\ k) 



6 Isospin matrix-elements calculation 



Independent of the functional form of the single-particle wave functions, the 
isospin dependence of the matrix elements is worked out in a proton-neutron 
basis representation. We find it convenient to introduce the proton/neutron 
creation/destruction operators, r±, in terms of its Cartesian components, as 



T± = 2 (^^ =t ^^J/) ' ^0 = 



(157) 



such that 



r+\p)= 
r_ \n) = 



r+|n)= \p) 
r_ \p) = \n) 
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Reciprocally, we have 



ry = - (r+ - r_ 



(158) 
(159) 



where the Cartesian components of the isospin operator are the usual Pauli 
matrices 



1 

1 



-t 

1 



1 
-1 



(160) 



With these definitions, the expectation values of the various isospin-dependent 
operators are computed, in terms of the matrices 



(r'|r+|r) 
{t I t_ I r) 

{r'\ro\T) 



1 if |r) = |n) and |r') = \p) , 

otherwise , 

1 if |r) = \p) and |r') = |n) , 
otherwise , 

1 if |r) = |r') = b) , 
-1 if |r) = |r') = \n) , 
otherwise . 



In particular, we note the identity: 

ri-T2 = 2 (ri+ra- + ri_r2+) + TioTso 



(161) 
(162) 

(163) 



(164) 



7 Results and Discussions 



To study the numerical accuracy of the approach described in the previous 
sections, we will consider all possible ph-coup\ed matrix elements, NME, of the 
Argonne fig potential that can be formed in a ^^0-like model space of single- 
particle harmonic-oscillator wave functions, T-LOni, with / < 6 and n < 6: the 
states involving nOio and nOu, i.e. (nlj) = {(10^), (11^), (ll|)}, are taken 
to be hole (occupied) states, and all others are particle (unoccupied) states. 
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Table 1 

Convergence of matrix elements of the central interaction of the Argonne vis po- 
tential calculated using the Moshinsky transformation-brackets approach outlined 
in Sec. [3l Here we study the evolution of the parameter e^^., which indicates how 
many matrix elements out of the total number of matrix elements, NME, that can 
be formed in a ^^0-like model space of single-particle wave functions, TZnh with / < 6 
and n < 6, are characterized by a relative error ranging between x% and y% when 
doubling the number of abscissas and weights in the Gauss-Hermite quadrature set, 
NGAUS, use to compute the radial integral ()32p . 

NGAUS NME e>i efl, ef,l, efZi 



64 23950 4792 16961 1907 270 

128 23950 122 8781 14595 

256 23950 64 7279 

512 23950 114 



Table 2 

Convergence of matrix elements of the spin-spin interaction of the Argonne vis 
potential calculated using the Moshinsky transformation-brackets approach outlined 



NGAUS 


NME 




<i 

^>0.1 


<o.i 

^>0.01 


<0.01 
^>0.001 


64 


43796 


14168 


20877 


8237 


476 


128 


43796 


172 


559 


16236 


24244 


256 


43796 





172 


725 


14238 


512 


43796 








179 


758 



We will begin by studying the convergence of the matrix elements calculated 
using the Moshinsky transformation-brackets approach outlined in Sec. |3l Up 
to intrinsic round-off errors, the numerical accuracy of the matrix element 
calculation is tied to the calculation of the radial part of the matrix element, 
given in Eq. ( I32l) . While this is not the only possible way to calculate this 
integral, we find it instructive to perform this integral using the Gauss-Hermite 
quadrature formula. This approach is applicable because the asymptotic form 
of the product of two harmonic-oscillator wave functions is Gaussian, and the 
potential itself falls to zero, as r goes to infinity. To study the convergence of 
the matrix elements with NGAUS, the number of abscissas and weights in the 
Gauss-Hermite quadrature set, it is convenient to introduce the parameter e^x, 
which represents the number of matrix elements characterized by a percentage 
change in the numerical value of the matrix element between x% and y% when 
doubling the number of grid points, NGAUS. 

Results for the 7 types of operators discussed in the Argonne fig potential are 
presented in Tables [IHZl We find that the calculation of the matrix-elements 
converges relatively fast with the number of points in the Gauss-Hermite 
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Table 3 

Convergence of matrix elements of the tensor interaction of the Argonne vis poten- 
tial calculated using the Moshinsky transformation-brackets approach outlined in 
Sec. [3l Similar to table [Tl 

NGAUS NME e>i efl, ef.j,, e|Hoi 



64 45853 16337 20967 7729 751 

128 45853 1907 12680 6663 3845 

256 45853 121 1851 12795 6667 

512 45853 6 147 2111 13314 



Table 4 

Convergence of matrix elements of the spin-orbit interaction of the Argonne vis 
potential calculated using the Moshinsky transformation-brackets approach outlined 
in Sec. El Simil ar to table [Tl 



NGAUS 


NME 




^>0.1 


<o.i 

^>0.01 


^<0.01 
^>0.001 


64 


41255 


7931 


20655 


10778 


1724 


128 


41255 


10 


108 


556 


5059 


256 


41255 





8 


78 


428 


512 


41255 








6 


66 



Table 5 

Convergence of matrix elements of the £^ interaction of the Argonne vis potential 
calculated using the Moshinsky transformation-brackets approach outlined in Sec.[3l 
Similar to table [Tl 



NGAUS 


NME 


e>i 


^>0.1 


<0.1 
^>0.01 


<0.01 
^>0.001 


64 


30650 


5560 


15907 


7636 


1421 


128 


30650 





130 


254 


7045 


256 


30650 








82 


158 


512 


30650 











82 



quadrature set, NGAUS, with the one notable exception of the matrix ele- 
ment of the tensor interaction. In order to converge the tensor matrix element 
we require a large set of Gauss-Hermite abscissas. Obtaining such a set is in 
fact a nontrivial endeavor, and we have described in |A] a practical approach 
to achieve this. The difficulty in converging numerically the tensor matrix el- 
ement is related to the fact that the radial shape of the tensor interaction 
has a longer range than the other components of the Argonne Wig potential. 
Hence the deviations from the tail of the integrand in Eq. (132|) are more 
pronounced in the case of the tensor interaction. A larger cutoff and a finer 
Gauss-Hermite grid is needed for an accurate numerical evaluation, which in 
turn leads to a larger NGAUS number of points. 
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Table 6 

Convergence of matrix elements of the £'^{ai ■ a2) interaction of the Argonne fig 
potential calculated using the Moshinsky transformation-brackets approach outlined 
in Sec. [3l Simil ar to table [T] 



NGAUS 


NME 




^>0.1 


<0.1 
^>0.01 


<0.01 
^>0.001 


64 


44261 


8280 


22408 


11709 


1715 


128 


44261 


2 


89 


765 


8686 


256 


44261 








2 


163 


512 


44261 











2 



Table 7 

Convergence of matrix elements of the quadrupole spin-orbit interaction in the 
Argonne vi^ potential calculated using Moshinsky transformation-brackets approach 
outlined in Sec . [3l Similar to table [T] 



NGAUS 


NME 




<i 

^>0.1 


<0.1 
^>0.01 


<0.01 
^>0.001 


64 


46726 


17572 


22013 


6245 


787 


128 


46726 


31 


408 


5549 


21721 


256 


46726 





10 


106 


1412 


512 


46726 








6 


49 



We conclude that a parameter value, NGAUS=512, is large enough to assure a 
good accuracy of all components of the Argonne wig matrix elements, including 
the tensor part of the interaction, as indicated by the sharp decrease in the 
parameter e>o.i for NGAUS=512. This statement is supported by results of 
an independent calculation, to be discussed next. 

Next, we will compare the calculation of the Argonne Wig matrix elements us- 
ing the Moshinsky transformation-brackets approach described in Sec. [3] with 
the more general calculation of the same matrix elements in the framework 
outlined in Sec. |5l With the exception of unavoidable round-off errors, the 
source of possible numerical loss of accuracy in the calculation of the matrix 
elements presented here is linked to the calculation of the integrals discussed 
in Sec. 15.11 First, we are concerned with the integrals in Eqs. f l75]) and flU^ : 
Given that the integrand vanishes for large values of the integration variables, 
and in order to maintain the approach general, we will evaluate these inte- 
grals by choosing a large cutoff and discretizing the integrand on a fine grid. 
Because these integrals only have to be performed once in the beginning of 
the calculation, and that the results can be stored on disk for future reference, 
we have chosen a radial cutoff of 25 fm and a momentum cutoff of 30 fm~^ in 
Eqs. fl73l) and fl92l) . respectively, together with uniform grids of 8000 points, 
and have performed these integrals using the Simpson rule. 

Second, we have to perform the integrals involving the radial parts of the 
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Table 8 

Convergence of matrix elements of the various components of the Argonne vis poten- 
tial calculated using the approach outlined in Sec.[5j Calculations were performed for 
values of the number of abscissas and weights in the Gauss-Hermite quadrature set, 
NGAUS, equal to 32, 64, 96 and 128. Here we study the evolution of the parameter 
e^^, which indicates how many matrix elements out of the total number of matrix 
elements, NME, that can be formed in a ^^0-like model space of harmonic-oscillator 
single-particle wave functions, TZni, with / < 6 and n < 6, are characterized by a rel- 
ative error between x% and y% when the number of grid points, NGAUS, changes. 
Calculations for larger values of NGAUS, (NGAUS=96 or 128 where appropriate) 
have resulted in zero values of the parameter e^^. 



Operator 


NGAUS 


NME 




^>0.1 


<0.1 
^>0.01 


<0.01 
^>0.001 


1 


64 


23950 





20 


411 


17243 


(Ji - (jj 


64 


43796 





137 


1425 


22048 


Sij 


64 


45853 


17844 


4704 


3795 


17892 


Sij 


96 


45853 





4 


23 


378 


i-S 


64 


41255 


10 


109 


636 


5191 


e 


64 


30650 





115 


257 


3011 




64 


44261 





9 


294 


3659 




64 


46726 


22 


230 


2971 


16467 



Table 9 

Comparison of the matrix elements of the 7 operators in the Argonne wis potential 
calculated using the center-of-mass separation of variables for harmonic-oscillator 
wave functions such as described in Sec. [3] and the general procedure outlined 
in Sec [5j Here, the parameter e^^ reflects the changes between the Moshinsky- 
transformation based result (NGAUS=512) and the matrix elements calculated us- 
ing the approach presented in Sec. [5] for NGAUS=96. 



Operator 


NME 


e>i 


<i 

^>0.1 


<0.1 
^>0.01 


<0.01 
^>0.001 


1 


23950 
















43796 











188 


Sij 


45853 





8 


187 


3365 


£-S 


41255 











2 


e 


30650 














{(y^■aj) 


44261 














{I ■ Sf 


46726 











4 



single-particle wave functions, as seen in Eqs. (!89|) and ( 19T]) . For an arbi- 
trary functional form of the single-particle wave functions, these calculations 
can become computationally expensive, as they may have to be repeated for 
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the calculation of each and every matrix element. For the case of a linear 
combination of harmonic-oscillator single-particle basis, these integrals can be 
performed using Gauss-Hermite quadrature formulas. 

We study the convergence of the matrix elements with the number of abscissas 
and weights in the Gauss-Hermite quadrature set, NGAUS, by computing the 
value of the parameter which represents the number of matrix elements 
characterized by a percentage change in the numerical value of the matrix 
element between x% and y% as the number of grid points, NGAUS, changes. 
We have performed calculations using 4 sets of Gauss-Hermite quadrature 
points, corresponding to values of the NGAUS parameter of 32, 64, 96, and 
128. Results for the 7 types of operators part of the Argonne vis potential 
are presented in Table [H We find that the calculation of the matrix-elements 
is fully converged for a value of the NGAUS parameter, NGAUS=96, for all 
components of the Argonne Vis potential, with the exception of the tensor in- 
teraction, which requires a larger set of NGAUS=128 Gauss-Hermite quadra- 
ture points. The slower convergence of the tensor-interaction matrix element 
follows the trend observed in the case of the calculation based on the Moshin- 
sky transformation brackets, and is due to the fact that the tensor interaction 
in the Argonne fig potential has a longer range and more structure than the 
other components of the interaction. 

In order to ascertain the numerical accuracy and correctness of the approach 
presented in Sec. we have also compared with results of the calculation 
of the matrix elements using the Moshinsky transformation brackets dis- 
cussed earlier. Results are displayed in Table [91 Here, the parameter e^^ 
reflects the changes between the Moshinsky transformation-brackets result 
(NGAUS=512) and the matrix elements calculated using the approach dis- 
cussed in Sec. \5\ for NGAUS=96. The two calculations lead to numerically 
identical matrix elements, which is quite remarkable since the two approaches 
are very different. This validates both theoretical approaches. 

In summary, in this paper we discussed two approaches for the calculation of 
two-body matrix elements of the Argonne fig potential. The first approach is 
only applicable to the specific case of a harmonic-oscillator single-particle wave 
function representation. In this case, the matrix elements are calculated us- 
ing the Talmi transformation (implemented numerically using the Moshinsky 
transformation brackets) which allows for the separation of the center-of-mass 
and relative coordinates degrees of freedom. Integrals involving the radial part 
of the potential were performed using Gauss-Hermite quadrature formulas, and 
convergence was achieved for sets of at least 512 Gauss points. This procedure 
was designed to validate the calculation of matrix elements of the Argonne fig 
potential using an approach suitable for an arbitrary functional form of the 
single-particle wave functions. The latter approach represents the main thrust 
of this paper. This general framework is suitable for the calculation of ma- 
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trix elements involving a representation of the single-particle wave functions 
given by linear combinations of harmonic-oscillator wave functions such as in 
Refs. pi?|llfl2] . and/or the two-scale functional representation for the hole 
and particle sides of the spectrum such as outlined in Ref. [2T]. For a model 
space represented in terms of harmonic-oscillator wave functions, results ob- 
tained using the approaches discussed in Sec. [3] and Sec. [5] are shown to be 
identical within numerical accuracy. 



A Generator of large- Gauss-Hermite quadrature sets 



The algorithm used to generate the Gauss-Hermite quadrature abscissas and 
weights represents an adaption of the algorithm described in Ref. 



The Gauss-Hermite quadrature formula allows for the efficient numerical com- 
putation of the integral 

/ dx fix) e-^ ^ ^3 fi^j) , (A.l) 

where Xj are the roots of the Hn{x) Hermite polynomial. The A^-points Gauss- 
Hermite becomes exact provided that f{x) is a polynomial. In order to avoid 
the overflow problems related to the standard recurrence relation of the Her- 
mite polynomials, i.e. 

Hj+,{x) = 2xHj{x) - 2jif,_i(x), (A.2) 

one uses instead the the orthonormal set of polynomials, Hj{x), which are 
generated using the recurrence relation 



#_i = 0, Ho = 7r-'/\ (A.3) 
Hj+iix) = X Hj{x) - Hj_i{x) . 

Then, the formula for the Gauss-Hermite weights is 

2 

Wi = , (A.4) 

where Hj{x) denotes the derivative of the polynomial Hj{x), which is calcu- 
lated as 

H;{x) = ^jH,.,ix). (A.5) 

The roots of the Hermite polynomial Hn{x) are symmetric about the origin, 
so we only have to calculate half of them, e.g. the positive ones. The roots are 
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estimated using the formulas put forward by Szego for the largest root of the 
Hermite polynomial Hjs[{x) [26], i.e. 



xj^' ~ V2iV + 1 - 1.85575 (2iV + (A.6) 
and by Stroud and Secrest [27] for the second-largest root, 

xn - l.UN-^'^yxN, (A.7) 



the third-largest root 

xj^La - 1.86X7V-1 - 0.86 xjv, (A.8) 

the fourth-largest root 

xS^Lg ^ 1.91 x^_2 - 0.91 (A.9) 

and all other positive roots 

xf ^ 2xj_i - Xj„2 • (A. 10) 

Subsequently, a root-finding algorithm is employed to improve these solutions. 
For added accuracy, all calculations are performed in quadruple precision, 
even though the Gauss- Hermite quadrature sets are ultimately used in double 
precision. 

Press et al. [25] propose the commonly-used Newton's method to achieve this 
task. In this approach, the i*^ approximation of the root Xj is obtained as 

^ = ^r^i - ^^^^44!^ . (A.ii) 

In our experience, this procedure is able to successfully produce all A^-points 
Gauss-Hermite quadrature sets of abscissas and weights for < 200. For 
larger sets, one must replace the root finder and use Muller's method instead 
(see page 364 in Ref. [25].) This method uses a quadratic interpolation among 
the 3 most recent estimates of the solution. To begin the search, we have 
chosen the first 3 guesses of the solution as xf \ given by the Szego's or Stroud's 

approximations described above, together with the coordinates xf^ ±0.1. It is 
important for this root-finding algorithm to be allowed to carry out the search 
in the complex plane, even though the roots of the Hermite polynomial are 
real. Muller's method allows for the calculation of Gauss-Hermite quadrature 
sets with A^ < 600. 

To further improve the algorithm and reach higher values of A^, it has been 
noticed that the coordinate dropped after the first estimate performed using 
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Muller's method, must be the one which is the farthest away from this esti- 
mate. This allows for the calculation of Gauss-Hermite quadrature sets with 
N < 1050, at which point we are very close to the natural numerical barrier 
encountered in numerical computations performed using 32-bit CPUs, even 
when using quadruple precision. Because, these sets where large enough to 
provide convergence for the purpose of our matrix-elements calculations, and 
because our calculations are carried out entirely on 32-bit platforms, we have 
not investigated this matter further. 

Because of the symmetry of the Gauss-Hermite quadrature set and the fact 
that our integrals involve only the radial coordinate in a spherical coordinates, 
we define the parameter NGAUS as half the value of the integer N in Eq. f lA.ip . 



B Reduced matrix elements 



In this appendix we list commonly used reduced matrix elements, relevant to 
the calculations described in this paper. 



B.l CW(f) 

The reduced matrix element of the unnormalized spherical harmonic is (|23] 
Eq. 5.4.6) 

{h II C('=)(f) II h) = 12 {hOkO I ZiO), (B.l) 
with the selection rules ii + £2 + k = even, and |/i — /2I < k < li + I2. 



B.2 \C^^^Xf)®V]^^^^ 



Using ([23], 7.1.1) we obtain 

{nM[C^'''\r)®^f'\\n2l2) (B.2) 
= (-)'^+'^+'^^2E r' ^ [(/i||C('=^)(^)IIO WI|V||n2/2) 

I [ ^2 ^1 ^ 

subject to the selection rules: |A;i — 1| < k2 < ki + 1. The reduced matrix 
elements of the gradient operator are obtained from the gradient formula (see 
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page 79 of Ref. [23]). The only nonzero components are 



(ni(/i = /2-l)||V||n2/2) 



and 



— roc 

Jo 



(B.3) 



= l2 + l)||V||n2/2) = V^2 + l_^" P dr] (r) 



d k 
dr r 



(B.4) 



B.3 [C(^i)(f) ® [V® V](^') 
Using ([23], 7.1.1) we obtain 



(fc2) 



(rii/i II \C^''\r) ® [V ® V](^')] ^"'^ 11^2/2) 

h j k2 



(B.5) 



'^+'^+^•^^2 ET" ' ^(^illC^'^H^)IIO WII[V®V](^')||n2/2) 

Z [ ^2 ^1 ^ 

with the selection rules: |£— 1| + 1. The allowed reduced matrix 

elements are 



= /2-2)||[V®V](^')||n2/2) = Jhik-l) 



1 1 3 
I2I2-2I2-I 



dr"^ r dr 



(B.6) 



Rn2h (r) ( , 



and 

(ni(/i = /2)||[V®V](^')||n2/2) 







> + {l2 + l)< 


fll . 1 




H 


11 













(B.7) 



x{-yj [r' dr] R„,Ur) (l)ir) 
Jo 



rf2 2d hik + 1) 
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and 



Mh=h + 2)\\[V^V]^^^\\n,k)= v'(/2 + l)(/2 + 2) \ ^ ^ ^ 

' I2I2 + 2I2 + I 



roo 

X {-Yj / pdr] (r) 

JO 



^2 2/2 + 1 /2(/2 + 2) 



(ir^ r dr 



(B.8) 
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